More Complex Mixed Models

Author

Marnin Wolfe

Previously: Intro to Mixed Models

Lesson plan

  1. Genetic Analysis with Mixed Models: introducing asreml package for more complex G and R structures compared with lmer. We’ll have a quick, non-comprehensive primer on asreml, followed by an example and a hands-on exercise.

ASReml Primer

Install asreml-R

First thing, you’ll need to install the asreml-R package. Please follow this written and video guide to download and install ASReml for R.

https://asreml.kb.vsni.co.uk/knowledge-base/asreml-r-installation-guide/

It’s not as simple as install.packages(). You’ll need the Activation Code I provided by email/Canvas. Note this is is a shared AU license including Horticulture, Forestry and CSES.

Introduction

Previously, we introduced mixed-models and compared them to linear models.

In lecture, we dug into the mechanics (matrix equations, variance-covariance structures).

There are lot’s of software options for mixed-models. They vary in their accessibility, ease-of-use, functionality and computational efficiency. It’s beyond our scope here, but an Ch. 6 of Ahmadi and Bartholomé (2022) gives a thorough overview as of 2022.

Here’s a great Introduction to ASReml-R (Tanaka 2018).

I highly recommend you reference the ASReml-R Reference Manual V4.2; it includes both code specifications AND theory. In my experience, to be successful, you’ll need to use it.

BE SURE YOU’RE ON V4.2 (V3 had very different syntax).

From Tanaka 2018 Slides

Estimates of fixed effects and prediction of random effects (Tanaka 2018)

The lmer() function we learned previously only handles the case when \(G\) and \(R\) matrices are \(I\). In other words when errors and random-effects are independent and identically distributed (IID). We’ll use asreml to access a richer universe of options.

The R user interface for ASReml (Tanaka 2018)

Basics

This section gives a quick primer on core asreml functionality. Examples with real data follow.

  • ASReml is a license-based propriety software.
  • It uses average information (AI)-REML to solve models.
  • The core software is based on FORTRAN and comes as a standalone program.
  • ASReml-R is an R-package that provides an interface to ASReml.
    • It is computationally pretty efficient thanks to using FORTRAN libraries under-the-hood rather than doing EVERYTHING in R.
    • One downside with using R: all datasets and matrices have to be held in your computers RAM. The smaller the RAM, the smaller the data you can handle.
asreml(respone ~fixed_formula,
       random = ~random_formula,
       residual = ~rcov_formula)

Variance Structures

For a complete list and specifications: Appendix C of the ASReml-R 4.2 Manual

Homogeneous Variances

# Example
fit1 <- asreml(Yield ~ Site,
         random= ~ Site:Geno,
         data=data1)

By default, asreml will assume random and residual terms are IID.

In the asreml manual, they call this the “homogeneous variance structure”.

Equivalently, you can write:

fit1a <- asreml(Yield ~ Site, data=data1
                random= ~ idv(Site):idv(Geno))

This stands for “identity variance” structure = \(var(u_{site})=\textbf{I}\sigma^2_{Site}\).

Same goes for the residual term. Here’s another equivalence:

fit1c <- asreml(Yield ~ Site, data=data1,
                random= ~ Site:Geno,
                residual= ~ idv(units))

NOTE: units is a special reserve word in asreml equivalent to factor(1:nrow(data1)). It only shows-up in residual terms.

Heterogeneous Variances

For various reasons, we might want to estimate a separate error-variance for each of several trials.

Or we might want to jointly analyze trials containing semi-overlapping populations / different numbers of genotypes. In such situations, you could postulate a different genetic variance for each site.

For this, we want “heterogeneous variance” structures.

In asreml there are a few ways to do this:

You can use idh(Site):Geno or diag(Site):Geno.

These fit separate variances for each level of the factor indicated within the ().

Conditional factors with at()

The at(f,l) specifics a structure that conditions level l and factor f. It means, we can choose to estimate a variance term only at certain levels of an effect. For example, if a one trial is replicated and another is not, you might want to do a joint analysis of both trials, including the replication term only for the replicated trial. E.g. at(Site,"PBU"):rep.

at() can be used in the G structure too. When the l (factor level) is not specified, e.g. at(Site), instead of at(Site,"PBU") the result is equivalent to diag() or idh(), a separate variance is estimated for each Site.

See ASReml-R manual 3.7.

Separable structures

Separability means that a covariance structure can be broekn down into the product of smaller matrices (non-technical definition, see ASReml-R manual 2.1.9).

\[ \Sigma_{u_{AB}} = var(u_{AB}) = \Sigma_A \otimes \Sigma_B \]

\(\otimes\) is called a “Kronecker” product in matrix algebra.

Illustration of a Kronecker Product from 10.48550/arXiv.2305.03875

It means that one variance structure can be “crossed” or “interacted” with another.

Here’s an example:

fit2 <- asreml(Yield ~ Site, 
              random=~ us(Site):id(Geno),
              data=data1)

The us(Site) stands for “unstructured” covariance. It means a difference variance and a different covariance are estimated for every level.

This one is useful for multi-trait models, among other things. WARNING: Convergence failures and matrix-inversion failures common as us() estimates all pairwise variance components. More terms to estimate = fewer degrees of freedom = higher amount of data needed.

becomes

Residual structures

In a multi-environment and multi-year trial, it is realistic to expect different amounts of experimental error across trials.

This situation can be handled using dsum(), which allows for Conditional Factors in the residual. Note that, in some cases, at() can be used equivalently in the residual.

dsum() “a direct sum of l variance matrices, one for each level of the conditioning factor… The data observations are partitioned into… sections to which separate variance structures are applied.” (see ASReml-R manual 3.7)

# generic
residual = ~dsum(~units | section)
# example
fit3 <- asreml(Yield ~ Site, 
              random=~ us(Site):id(Geno),
              residual = ~dsum(~units|Site),
              data=data1)

dsum() can also specify different structures on the residuals within different sections (e.g. different Sites.

RULES for Residual Term:

  1. Number of effects in residual must be equal to number of observations.
  2. Compound terms in residual need to ensure each combo of levels uniquely identifies one unit of data.
  3. The data must be ordered to match the residual. E.g. sort observations within Sites for ~dsum(units | Site).

asreml Output

Basic asreml output:

# Variance Parameters
summary(fit1)$varcomp # or use lucid::vc(fit1)

##  BLUPs / Prediction

# E-BLUE
coef(fit1)$fixed

# E-BLUP
coef(fit1)$random

LRTs and model comparison

Recall the likelihood ratio test discussed previously:

\[ LRT = -2 \times ln \left(\dfrac{L(reduced)}{L(full)}\right) \]

This ratio is compared to a \(\chi^2\)-distribution with degrees of freedom equal to the diff. in number of variance components between the two models.

A small p-value for the test indicates the two models are different.

It DOES NOT necessarily mean your “full” model is the better model. Check the likelihood, or compare the AIC/BIC. If the full model has higher likelihood, then declare the focal variance component “significant”.

Information Criterion

Akaike Information Criterion (AIC), and the Bayesian Information Criterion (BIC) are related measures of the goodness-of-fit for a model.

\[AIC=-2 \times logL + 2 \times p\]

\[BIC = -2 \times logL + 2 \times p \times log(n)\] - p is the number of varaince parameters in the model - n the number of observations

With AIC and BIC, lower-is-better, i.e. the smaller AIC is the better model.

Both of these criteria prioritize models with a high likelihood while penalizing model complexity.

# Fixed effects 
# Wald test
fit <- asreml(yield ~ Variety*Nitrogen, random= ~ Blocks/Wplots, data=oats)

wald(fit)

# Rand effects
full_model <- asreml(yield ~ gen, random=~rep, trace=F,
               data=burgueno.rowcol)
null_model <- asreml(yield ~ gen, trace=F,
               data=burgueno.rowcol)
lrt(fit2, fit1)

Other

Convergence issues

  • Run more iterations:
    • increase maxiter=
    • Use update(fit)
    • Use initial values (init) in the argument of the asreml function

Memory space

  • Increase the RAM available to ASREML (workspace="10gb")
  • Use top, htop, Activity/System Monitor –> watch you system resources

Example

Let’s keep using the data subset we chose for previous lessons.

Data

library(tidyverse)
library(asreml)
Online License checked out Wed Jan 28 11:01:48 2026
library(lme4)
alt<-read.csv(here::here("data","19-25CC_ALT_DATA_ALL_2025-08-23.csv"))

# My example data chunk, same as before
# Keep using your own!
al_alt<-alt %>% 
  filter(site=="AL")

# Make factors before modeling
al_alt<-al_alt %>% 
  mutate(entry=as.factor(entry),
         site.year=as.factor(site.year),
         # explicitly nest values of rep in site.year
         # review explicit vs. implicit nesting
         rep=as.factor(paste0(site.year,"-",rep)))

lmer vs. asreml

The lmer model we fit was:

lmer_mm<-lmer(biomass.1 ~site.year + (1|rep) + (1|entry) + (1|entry:site.year),
         data = al_alt)
summary(lmer_mm)
Linear mixed model fit by REML ['lmerMod']
Formula: 
biomass.1 ~ site.year + (1 | rep) + (1 | entry) + (1 | entry:site.year)
   Data: al_alt

REML criterion at convergence: 939.3

Scaled residuals: 
     Min       1Q   Median       3Q      Max 
-2.18776 -0.60236  0.00158  0.57017  2.76647 

Random effects:
 Groups          Name        Variance Std.Dev.
 entry:site.year (Intercept)  58.935   7.677  
 entry           (Intercept)  17.335   4.163  
 rep             (Intercept)   5.807   2.410  
 Residual                    142.962  11.957  
Number of obs: 116, groups:  entry:site.year, 58; entry, 35; rep, 4

Fixed effects:
              Estimate Std. Error t value
(Intercept)     50.027      2.858  17.506
site.year25AL  -11.592      3.878  -2.989

Correlation of Fixed Effects:
            (Intr)
site.yr25AL -0.694

We can fit the equivalent model with asreml like so:

asr_mm<-asreml(biomass.1 ~site.year,
               random = ~entry + rep + site.year:entry,
               data = al_alt)
ASReml Version 4.2 28/01/2026 11:01:49
          LogLik        Sigma2     DF     wall
 1     -366.2039      175.6047    114   11:01:49
 2     -365.4699      164.5684    114   11:01:49
 3     -364.8944      145.4173    114   11:01:49
 4     -364.8831      143.0244    114   11:01:49
 5     -364.8830      142.9635    114   11:01:49
summary(asr_mm)$varcomp
                 component std.error   z.ratio bound %ch
rep               5.806474  10.86064 0.5346348     P 0.1
entry            17.335959  32.80183 0.5285059     P 0.1
site.year:entry  58.934355  40.82486 1.4435898     P 0.0
units!R         142.963533  27.02516 5.2900159     P 0.0
# lucid::vc(asr_mm) ## this just makes a prettier looking output
# E-BLUE
coef(asr_mm)$fixed
                  effect
(Intercept)     50.02667
site.year_24AL   0.00000
site.year_25AL -11.59152
# E-BLUP for entry
coef(asr_mm,list = T)$entry %>% head
                       effect
entry_17MDCCH       1.0255866
entry_17MDCCS       2.5387950
entry_17RMCCL      -1.8767066
entry_18NCCCEGiant  1.3787976
entry_19MDCC        1.8343968
entry_19NCCCLH     -0.3611015

This gives us just the BLUPs.

# Also E-BLUP
summary(asr_mm,coef=T)$coef.random %>% head
                     solution std.error     z.ratio
entry_17MDCCH       1.0255866  3.718116  0.27583502
entry_17MDCCS       2.5387950  3.718116  0.68281758
entry_17RMCCL      -1.8767066  3.922077 -0.47849816
entry_18NCCCEGiant  1.3787976  3.921467  0.35160252
entry_19MDCC        1.8343968  3.718116  0.49336729
entry_19NCCCLH     -0.3611015  3.718116 -0.09711949

This gives us also the standard errors, which when squared give us the Prediction Error Variances (PEVs).

Check assumptions

We can plot residuals and check assumptions in much the same was as with plot.lmer.

plot(asr_mm)
Warning: Using `size` aesthetic for lines was deprecated in ggplot2 3.4.0.
ℹ Please use `linewidth` instead.
ℹ The deprecated feature was likely used in the asreml package.
  Please report the issue to the authors.

Hypothesis testing with asreml

Now let’s do a likelihood ratio test.

asr_nullentry<-asreml(biomass.1 ~site.year,
                      random = ~rep + site.year:entry,
                      data = al_alt)
ASReml Version 4.2 28/01/2026 11:01:49
          LogLik        Sigma2     DF     wall
 1     -367.3448      192.4790    114   11:01:49
 2     -366.0755      174.8798    114   11:01:49
 3     -365.2521      156.9297    114   11:01:49
 4     -365.0124      143.9875    114   11:01:49
 5     -365.0110      142.9698    114   11:01:49
lrt(asr_mm,asr_nullentry)
Likelihood ratio test(s) assuming nested random models.
(See Self & Liang, 1987)

                     df LR-statistic Pr(Chisq)
asr_mm/asr_nullentry  1      0.25599    0.3064

That’s a bummer. The genetic effect in the data chunk does not appear significant.

Test on the fixed effects

wald(asr_mm)
Wald tests for fixed effects.
Response: biomass.1
Terms added sequentially; adjusted for those above.

              Df Sum of Sq Wald statistic Pr(Chisq)    
(Intercept)    1     65703         459.58 < 2.2e-16 ***
site.year      1      1278           8.94  0.002792 ** 
residual (MS)          143                             
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

NOTE: if you really want to do sig. tests on your mixed-model’s fixed effects, do some reading first. There are some options we should consider regarding the Wald test. Namely, I believe it does an equivalent of a Type I SS (sequential). I think (am not certain) you can achieve something like a Type III (non-sequential) with something like this wald(fit, ssType="conditional", denDF="numeric"). Verify.

Variance structures

Next, let’s fit a more appropriate data to the model, one that lmer cannot do.

There are several possibilities: 1. Fit a heterogeneous (different) error variance for each site.year 2. There are a few entry in 24AL that are not in 25AL, one could estimate a different genetic variance for each site.year. 3. Variance among complete blocks could differ between site.years.

For demo purposes, I’ll focus on the heterogeneous errors. That’s what’s most plausible IMHO.

asr_heterror<-asreml(biomass.1 ~site.year,
                     random = ~entry + rep + site.year:entry,
                     residual = ~dsum(~units|site.year),
                     data = al_alt %>% arrange(site.year,X))
ASReml Version 4.2 28/01/2026 11:01:49
          LogLik        Sigma2     DF     wall
 1     -371.5225           1.0    114   11:01:49
 2     -367.8455           1.0    114   11:01:49
 3     -364.3065           1.0    114   11:01:49
 4     -362.4049           1.0    114   11:01:49
 5     -362.2450           1.0    114   11:01:49
 6     -362.2319           1.0    114   11:01:49
 7     -362.2270           1.0    114   11:01:49
 8     -362.2242           1.0    114   11:01:49
 9     -362.2221           1.0    114   11:01:49
10     -362.2202           1.0    114   11:01:49  (  1 restrained)
11     -362.2197           1.0    114   11:01:49  (  1 restrained)
12     -362.2197           1.0    114   11:01:49  (  1 restrained)
# summary(asr_heterror)$varcomp
lucid::vc(asr_heterror)
           effect component std.error z.ratio bound %ch
              rep   0.00001        NA      NA     B  NA
            entry  20.7         37.16    0.56     P   0
  site.year:entry  74.16        43.06    1.7      P   0
 site.year_24AL!R 196           47.75    4.1      P   0
 site.year_25AL!R  80.66        21.29    3.8      P   0

Notice that two separate error variance estimates are produced, one for each site.year.

Looks like 25AL had smaller error than 24AL.

Now we have a nested model between the original and the one with added het. errors.

Two ways to compare:

lrt(asr_mm,asr_heterror)
Likelihood ratio test(s) assuming nested random models.
(See Self & Liang, 1987)

                    df LR-statistic Pr(Chisq)  
asr_heterror/asr_mm  1       5.3265    0.0105 *
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Ok, so the LRT says these two models are different.

plot(asr_heterror)

With complex residual structures, it appears asreml won’t auto-plot your residual diagnostics.

Distribution of residuals

hist(residuals(asr_heterror))

Standardized residuals vs. fitted.

plot(x=fitted(asr_heterror),y=scale(residuals(asr_heterror))); abline(b = 0,h = 0)

QQ plot:

res<-residuals(asr_heterror)
qqnorm(res); qqline(res, col = "red")

Model comparison

AIC and BIC are computed by asreml and available in the summary object. As long as you’ve got the same data and fixed-effects, you can use this along with the LRT to determine which model is best fit.

summary(asr_mm)$aic
[1] 737.766
attr(,"parameters")
[1] 4
summary(asr_heterror)$aic
[1] 734.4394
attr(,"parameters")
[1] 5

So the heterror model does have a lower AIC value. Along with the LRT p-value, we can declare this a significantly better model.

asr_heterror_nullg<-asreml(biomass.1 ~site.year,
                           random = ~rep + site.year:entry,
                           residual = ~dsum(~units|site.year),
                           data = al_alt %>% arrange(site.year,X), 
                           maxit=90)
ASReml Version 4.2 28/01/2026 11:01:50
          LogLik        Sigma2     DF     wall
 1     -376.5979           1.0    114   11:01:50
 2     -371.1740           1.0    114   11:01:50
 3     -365.8315           1.0    114   11:01:50
 4     -363.1212           1.0    114   11:01:50
 5     -362.4069           1.0    114   11:01:50
 6     -362.3648           1.0    114   11:01:50
 7     -362.3575           1.0    114   11:01:50
 8     -362.3539           1.0    114   11:01:50
 9     -362.3515           1.0    114   11:01:50
10     -362.3495           1.0    114   11:01:50  (  1 restrained)
11     -362.3482           1.0    114   11:01:50  (  1 restrained)
12     -362.3481           1.0    114   11:01:50  (  1 restrained)
13     -362.3481           1.0    114   11:01:50  (  1 restrained)
lrt(asr_heterror_nullg,asr_heterror)
Likelihood ratio test(s) assuming nested random models.
(See Self & Liang, 1987)

                                df LR-statistic Pr(Chisq)
asr_heterror/asr_heterror_nullg  1       0.2568    0.3062

Yeah, well, still not showing significant genetic variation. Bummer.

summary(asr_heterror)$aic - summary(asr_heterror_nullg)$aic
[1] 1.743201
attr(,"parameters")
[1] 5

The AIC of the full model is higher than the null.

What about the GxE term?

asr_heterror_nullgxe<-asreml(biomass.1 ~site.year,
                           random = ~entry + rep,
                           residual = ~dsum(~units|site.year),
                           data = al_alt %>% arrange(site.year,X), 
                           maxit=90)
ASReml Version 4.2 28/01/2026 11:01:50
          LogLik        Sigma2     DF     wall
 1     -376.8922           1.0    114   11:01:50
 2     -372.0801           1.0    114   11:01:50
 3     -367.5884           1.0    114   11:01:50
 4     -365.4995           1.0    114   11:01:50
 5     -364.8730           1.0    114   11:01:50
 6     -364.7547           1.0    114   11:01:50  (  1 restrained)
 7     -364.7217           1.0    114   11:01:50  (  1 restrained)
 8     -364.7161           1.0    114   11:01:50  (  1 restrained)
 9     -364.7154           1.0    114   11:01:50  (  1 restrained)
lrt(asr_heterror_nullgxe,asr_heterror)
Likelihood ratio test(s) assuming nested random models.
(See Self & Liang, 1987)

                                  df LR-statistic Pr(Chisq)  
asr_heterror/asr_heterror_nullgxe  1       4.9914   0.01274 *
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
summary(asr_heterror)$aic - summary(asr_heterror_nullgxe)$aic
[1] -2.991408
attr(,"parameters")
[1] 5

Ok! So the GxSiteYear term is significant and improves model fit.

There’s more that can be done here, alternative models that could be tested.

Good enough for our demo.

PEV and REL

Let’s say we are happy with our model.

Now we want to extract BLUPs, and compute the Reliability of Prediction.

The summary() function includes the “standard errors” of prediction.

blups<-summary(asr_heterror,coef=T)$coef.random %>% 
  # confert to data.frame
  as.data.frame %>% 
  # change the rownames to a column the the data.frame
  rownames_to_column(var = "entry") %>% 
  # filter to only the BLUPs for the entry main-effect 
  # (output includes all random terms)
  filter(grepl("^entry_",entry)) %>% 
  # for good measure, tidy data, 
  # need to remove a prefix that asreml added to each entry name
  mutate(entry=gsub("entry_","",entry,fixed = T))
blups %>% dim
[1] 35  4
blups %>% head
         entry   solution std.error    z.ratio
1      17MDCCH  1.5784997  4.008503  0.3937878
2      17MDCCS  2.7495080  4.008503  0.6859189
3      17RMCCL -1.7051249  4.307455 -0.3958544
4 18NCCCEGiant  1.7966818  4.199167  0.4278663
5       19MDCC  2.1806776  4.008503  0.5440129
6     19NCCCLH -0.5985299  4.008503 -0.1493151

Recall that \[REL_i = 1 - \frac{PEV_i}{\sigma^2_u}\]

Vg<-summary(asr_heterror)$varcomp["entry","component"]
blups<-blups %>% 
  mutate(PEV=std.error^2,
         REL=1-(PEV/Vg))

hist(blups$REL)

Recall that the average reliability of prediction for a main genetic effect is an estimator of the heritability (ratio of genetic to total phenotypic variance) (Cullis, Smith, and Coombes (2006); Schmidt et al. (2019)) .

mean(blups$REL)
[1] 0.1914108

Heritability

We haven’t talked much about Quantitative Genetics parameters yet.

The heritability is a critical parameter in genetics. It expresses that fraction of trait variation that is genetic. It is used, among other things, to predict the response of a population to selection. See “Breeder’s Equation”; we’ll talk about it soon!

Let’s compare the Cullis estimate (~0.19) to the “traditional” variance component based one.

The formula should be something like this:

\[ H^2 = \frac{\sigma^2_g}{\sigma^2_g + \sigma^2_{g\times yr} + \frac{\bar{\sigma^2_{e_k}}}{nYrs}} \]

Note I took the mean error across site.years.

varcomps<-summary(asr_heterror)$varcomp
VgXsiteyr<-varcomps["site.year:entry","component"]
Ve_meanPerSite<-mean(varcomps[grepl("site.year_",rownames(varcomps)),"component"])

Vg / (Vg + VgXsiteyr + Ve_meanPerSite)
[1] 0.08875425

Hands-on

Now it’s your turn!

  • Our goal is to be able to jointly analyze the entire ALT dataset and to produce estimates of genetic parameters (heritability) and rankings of the entries (BLUPs).
  • OPTIONS:
    1. Use your previous data chunk (important it includes multiple site-years)
    2. Pick a bigger, more complex chunk
    3. For the very bold: try analyzing the entire dataset.
  • Postulate and write the model you will use.
  • Find the best-fit model for the dataset
    • Determine this using Likelihood Ratio Tests and the AIC/BIC
  • Is there “significant” genetic variability in your data? Which variance components are >0?
  • What is the heritability?
    • Compute both the classical and the Cullis (mean REL) heritability.
  • Other things that should be instructive:
    • Compare you BLUPs to BLUEs and/or to raw means.
    • Explore the relationship between amount-of-data (number of plots, years, etc) and the REL of prediction. Why do some entries have lesser or greater reliability?
  • Which is the “best” Crimson Clover entry in the bunch? Does it vary by year? Location? Region?

Next

Advanced Mixed Model Topics: Probable Topics (multi-env. trial analysis, two-stage analysis, GxE, spatial models, heritability)

Hands-on Data Wrangling and Exploratory Data Analysis: how to manipulate and plot your data, functions, loops and analyzing multiple traits.

References

Ahmadi, Nourollah, and Jérôme Bartholomé, eds. 2022. Genomic Prediction of Complex Traits. Springer US. https://doi.org/10.1007/978-1-0716-2205-6.
Cullis, B. R., A. B. Smith, and N. E. Coombes. 2006. “On the Design of Early Generation Variety Trials with Correlated Data.” Journal of Agricultural, Biological, and Environmental Statistics 11 (4): 381–93. https://doi.org/10.1198/108571106x154443.
Schmidt, Paul, Jens Hartung, Jörn Bennewitz, and Hans-Peter Piepho. 2019. “Heritability in Plant Breeding on a Genotype-Difference Basis.” Genetics 212 (4): 991–1008. https://doi.org/10.1534/genetics.119.302134.